R Notebook

This is for the latest 2019 PLANEA test!

Easy ggplot theme.

Let’s define some function to interact easily with the database.

#' @title plan_connect
#'
#' @description This function takes the Postgres environment variables
#' to automatically connect to the planeadb database.
#' The connection must be saved into a variable in order to
#' use future functions. I'm cheating and not reading from .env.
#'
#' @examples con <- prev_connect()
#'
#' @param No parameter is needed.
#' @export
plan_connect <- function(){
  DBI::dbConnect(RPostgreSQL::PostgreSQL(),
  host     =  "localhost",
  user     =  "planea",
  password =  "planea",
  port     =  5432,
  dbname   =  "planea")
}

#' @title load_query
#'
#' @description Allows to run specific queries with the dbrsocial syntax.
#'
#' @param connection DBI connection. A connection to a database
#' @param schema variable. A valid schema from a database on the
#' @param the_table.  An existing table in the given schema.
#' @param colums string. The columns in the database we want to retrieve
#' information
#' @param options string. Part of the SQL query with containing WHERE, ORDER,
#' LIMIT and so statements
#'
#' @examples geom_muni <-
#' load_query(con,raw,sifode,columns="entidadfederativa",options="WHERE
#' mes='Abril'")
#' @export
load_query <- function(connection,schema,the_table,columns="*",options=""){
    the_query <- "SELECT %s FROM %s.%s"
    complete <- paste0(the_query," ",options)
    schema    <- deparse(substitute(schema))
    the_table <- deparse(substitute(the_table))
    initial <- RPostgreSQL::dbSendQuery(connection,
                             sprintf(complete,columns,schema,the_table))
    return(initial)
}

#' @title load_table
#'
#' @description This function loads a connection to
#' a given table from a particular schema in a database
#' connection.
#'
#' @param connection DBI connection. A connection to a database
#' must be open and given.
#' @param schema variable. A valid schema from a database on the
#' connected database.
#' @param the_table. An existing table in the given schema.
#'
#' @examples the_dic<-load_table(con,raw,sifode_dic)
#' @export
load_table <- function(connection,schema,the_table){
    the_query <- "SELECT * FROM %s.%s"
    schema    <- deparse(substitute(schema))
    the_table <- deparse(substitute(the_table))
    initial <- RPostgreSQL::dbSendQuery(connection,
                             sprintf(the_query,schema,the_table))
}

#' @title large_table
#'
#' @description This function loads a connection to a large table without
#' loading it to memory.
#'
#' @param connection DBI connection. A connection to a database
#' must be open and given.
#' @param schema variable. A valid schema from a database on the
#' connected database.
#' @param the_table. An existing table in the given schema.
#'
#' @examples cuis_table <- large_table(con,raw,cuis_39_9)
#' @export
large_table <- function(connection,schema,the_table){
    schema    <- deparse(substitute(schema))
    the_table <- deparse(substitute(the_table))
    retrieved <- dplyr::tbl(connection,dbplyr::in_schema(schema,the_table))
    return(retrieved)
}

#' @title clear_results
#'
#' @description This function clears the results from a previous executed
#' query.
#'
#' @param connection DBI connection. A connection to a database must be open and given.
#'
#' @examples clear_results(con)
#' @export
clear_results <- function(connection){
    DBI::dbClearResult(DBI::dbListResults(connection)[[1]])
}

#' @title join_tables
#'
#' @description Returns a match between two tbl by defined key
#' @param left_table a tbl-like object
#' @param right_table a tbl-like object
#' @param left_key the column name from left_table to compare
#' @param right_key the column name from right_table to compare
#'
#' @examples join_tables(cuis_sample,llave_hogar_h,domicilios_sample_query,llave_hogar_h)
#' @export
join_tables <- function(left_table, left_key, right_table, right_key){
    left_key <- substitute(left_key)
    right_key <- deparse(substitute(right_key))
    in_tables <- left_table %>%
        dplyr::filter(left_key %in% right_table[[right_key]])
    return(in_tables)
}

#' @title retrieve_result
#'
#' @description Return the fetch results of a query
#' @param query An exec unfetched query
#'
#' @examples sample_table(load_table(prev_connect(),raw,sifode))
#' @export
retrieve_result <- function(query,n=-1,number=Inf){
    if (class(query)[1] == "tbl_dbi"){
        the_table <- dplyr::collect(query,n=number)
        return(the_table)
    }
    else{
    the_table <- DBI::dbFetch(query,n)
    return(the_table)
    }
}

#' @title load_geom
#'
#' @description Gives a "ready to go" data frame for geometry plotting
#'
#' @param connection DBI connection. A connection to a database
#' @param schema variable. A valid schema from a database on the
#' @param the_table.  An existing table in the given schema.
#' @param colums string. The columns in the database we want to retrieve.
#' Defaults cve_ent, cve_mun, cve_muni.
#' @param geom_col . The name of the column in the database that contains a geometry
#' @param col_shape. The name of the column that we want to use to join
#' information
#' @param options string. Part of the SQL query with containing WHERE, ORDER,
#' LIMIT and so statements
#'
#' @examples geom_muni <- load_geom(con1,raw,geom_municipios,geom_col=geom,col_shape=cve_muni,options=options)
#' @export
load_geom <- function(connection,schema,the_table,columns="\"CVEGEO\", \"CVE_ENT\", \"CVE_MUN\"", geom_col, col_shape, options=""){
    geom_col <- deparse(substitute(geom_col))
    schema    <- deparse(substitute(schema))
    the_table <- deparse(substitute(the_table))
    col_shape <- deparse(substitute(col_shape))
    geom_col2 <- paste0("\"",geom_col,"\"")

    the_query <- "SELECT %s FROM %s.%s"
    geom_col_as <- sprintf(", %s as geom",geom_col2)
    columns <- paste0(columns,geom_col_as)
    complete <- paste0(the_query," ",options)

    initial <- RPostgreSQL::dbSendQuery(connection,
                             sprintf(complete,columns,schema,the_table)) %>%
    retrieve_result()

    mun_shp = rangeMapper::WKT2SpatialPolygonsDataFrame(initial, geom = "geom", id = col_shape)
    mun_df <- fortify(mun_shp, region = col_shape)
    names(mun_df)[names(mun_df)=="id"] <- col_shape

    return(mun_df)
}

#' @title discon_db
#'
#' @description This function disconnects a PostgreSQL
#' connection.
#'
#' @param connection DBI connection. A connection to a database must be open and given.
#'
#' @examples discon_db(con)
#' @export
discon_db <- function(connection){
    RPostgreSQL::dbDisconnect(connection)
}

First we load the PLANEA results.

We load the geometries but will use them later

planea[logros] <- sapply(planea[logros],as.numeric)

NAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercion

We need to change some municipalities names in order to join this database with the CONAPO one. This are the missed-matched names in the tables of the database. I’ll change the ones from the schools.

planea$municipio <- plyr::mapvalues(planea$municipio, to=to_conapo, from=from_schools)

The following `from` values were not present in `x`: DR. BESLISARIO DOMINGUEZ

Load the CONAPO table

Load the marginalization table and change the variable names to usable ones.

clear_results(con)

[1] TRUE

Join the three tables. I should run this in Postgres. TODO write SQL statement.

This result is for a table in the document

How many schools

The number of schools in a municipality.

conapo$CLAVE %>% unique() %>% length()

[1] 2457

Fit porwer-law to number of schools

schools_power_law <- lm(log10(num_school) ~ log10(number),num_schools)

schools_power_law

Call:
lm(formula = log10(num_school) ~ log10(number), data = num_schools)

Coefficients:
  (Intercept)  log10(number)  
        3.007         -1.389  

the_p_law <- tibble(number=seq_schools,value=adjusted_power)
num_schools %<>% left_join(the_p_law,by="number")

Number of evaluated students in language

sum(as.numeric(planea$evaluados_len))

[1] 2032174

Number of evaluated students in mathematics

sum(as.numeric(planea$evaluados_mat))

[1] 2024534

There’s a variable called representativity for both math and language but it’s unclear the percentage cut-off, then I will use the schools with above 50% of both tests

schools <- schools[(schools$porc_len > 50 & schools$porc_mat > 50),]
mun_schools <- schools %>%
  select(CLAVE,escuela,POP,marginacion,gm,elementary,illiteracy) %>%
  unique() %>% dplyr::group_by(CLAVE) %>%
  dplyr::summarise(num_school = n(),POB=unique(POP),elementary=unique(elementary),illiteracy=unique(illiteracy),gm=unique(gm)) %>%
  mutate(spc = num_school/(POB/100000))

Number of evaluated students in language

sum(as.numeric(schools$evaluados_len))

[1] 4032596

Number of evaluated students in mathematics

sum(as.numeric(schools$evaluados_mat))

[1] 4018728

unique((schools$clave)) %>% length()

unique((schools$CLAVE)) %>% length()

schools_l <- schools %>% group_by(CLAVE, clave, escuela, gm) %>% summarise(num_school=n(),POP=unique(POP),illiteracy=unique(illiteracy),elementary=unique(elementary),I_len = median(I_porc_len),II_len = median(II_porc_len),III_len = median(III_porc_len),IV_len = median(IV_porc_len),
                                                    I_mat = median(I_porc_mat),II_mat = median(II_porc_mat),III_mat = median(III_porc_mat),IV_mat = median(IV_porc_mat))

`summarise()` has grouped output by 'CLAVE', 'clave', 'escuela'. You can override using the `.groups` argument.

Translate some variables to english. There are some with peculiar name that shouldn’t be translated.

Define the colors

altoc <- "#e3b23c"
bajoc <- "#725752"
medioc <- "#fe64a3"
maltoc <- "#78bc61"
mbajoc <- "#4f359b"

Define plot colors

comc <- "#efecca"
privc <- "#db2b39"
pubgc <- "#f7ff58"
pubtc <- "#ff934f"
telec <- "#29335c"

schools with 100% in IV level for mathematics

Load the indigenous population and joint information with the test.

Plot the percentage of indigenous population by level of marginalization.

Some variables that won’t make it into the document

As expected the ratio of indigenous population in a municipality correlates perfectly with the percentage of inhabitants that speak a indigenous language.

The more indigenous more illiteracy.

The percentage of indigenous population correlates possitevily with the percentage of level 1 achievement.

lm(I_len/100 ~ p_len, mar_scores)

While for the level IV it has a negative correlation

lm(IV_len/100 ~ p_len, mar_scores)

How does language and math level I correlate in a highly marginalized municipality?

what about a very low marginalized area?

lm((I_mat) ~ (I_len),mar_scores[mar_scores$gm == "Very Low",])

Call:
lm(formula = (I_mat) ~ (I_len), data = mar_scores[mar_scores$gm == 
    "Very Low", ])

Coefficients:
(Intercept)        I_len  
    29.0414       0.8352  

lm((I_mat) ~ (I_len),mar_scores[mar_scores$gm == "Very High",])

Call:
lm(formula = (I_mat) ~ (I_len), data = mar_scores[mar_scores$gm == 
    "Very High", ])

Coefficients:
(Intercept)        I_len  
    32.9009       0.6141  

mt <- ggplot(mar_scores)+
  geom_point(aes(log10(IV_len),log10(IV_mat),color=gm))+
  scale_color_manual(name = "Marginalization", values = c("High" = altoc, "Low" = bajoc, "Medium" = medioc, "Very High" = maltoc, "Very Low" = mbajoc),
                      labels = c("Very Low", "Low", "Medium", "High", "Very High"))+
  coord_fixed()+
  geom_abline()+
  scale_x_log10(name = "percentage language")+#,limits = c(0,100))+
  scale_y_log10(name = "percentage mathematics")+#,limits = c(0,100))+
  ggtitle("Achievement Level IV by marginalization")+
  th+
  theme(
    axis.text.x = element_text(angle = 0),
    panel.grid.minor = element_blank(),
    legend.position = "none")

mt + facet_grid(cols=vars(gm)) +
    theme(strip.background = element_rect(fill="#b53a31"))

ggsave("../../figs/math_vs_language_IV_log.png",height = 20, width = 40, units = "cm")

mt <- ggplot(mar_scores)+
  geom_point(aes((IV_len),(IV_mat),color=gm))+
  scale_color_manual(name = "Marginalization", values = c("High" = altoc, "Low" = bajoc, "Medium" = medioc, "Very High" = maltoc, "Very Low" = mbajoc),
                      labels = c("Very Low", "Low", "Medium", "High", "Very High"))+
  coord_fixed()+
  geom_abline()+
  scale_x_continuous(name = "percentage language")+#,limits = c(0,100))+
  scale_y_continuous(name = "percentage mathematics")+#,limits = c(0,100))+
  ggtitle("Achievement Level IV by marginalization")+
  th+
  theme(
    axis.text.x = element_text(angle = 0),
    panel.grid.minor = element_blank(),
    legend.position = "none")

mt + facet_grid(cols=vars(gm)) +
    theme(strip.background = element_rect(fill="#b53a31"))

lm((IV_mat) ~ (IV_len),mar_scores[mar_scores$gm == "Very Low",])
lm((IV_mat) ~ (IV_len),mar_scores[mar_scores$gm == "Very High",])

mt <- ggplot(mar_scores)+
  geom_point(aes((IV_len),(I_mat),color=gm))+
  scale_color_manual(name = "Marginalization", values = c("High" = altoc, "Low" = bajoc, "Medium" = medioc, "Very High" = maltoc, "Very Low" = mbajoc),
                      labels = c("Very Low", "Low", "Medium", "High", "Very High"))+
  coord_fixed()+
  geom_abline(intercept = 100, slope=-1)+
  # scale_x_continuous(name = "percentage language",limits = c(0,100))+
  # scale_y_continuous(name = "percentage mathematics",limits = c(0,100))+
  ggtitle("Achievement Level IV by marginalization")+
  th+
  theme(
    axis.text.x = element_text(angle = 0),
    panel.grid.minor = element_blank(),
    legend.position = "none")

mt + facet_grid(cols=vars(gm)) +
    theme(strip.background = element_rect(fill="#b53a31"))

# ggsave("~/Documents/UVM/courses/DS2/assignments/HW6/math_vs_language_IV.png",height = 20, width = 40, units = "cm")

mt <- ggplot(mar_scores)+
  geom_point(aes((prop_indi),(IV_len),color=gm),alpha=0.6)+
  # geom_density(aes(pm_analf),)+
  # geom_point(aes(ph_analf,prop_indi,color=gm),alpha=0.4)+
  scale_color_manual(name = "Marginalization", values = c("High" = altoc, "Low" = bajoc, "Medium" = medioc, "Very High" = maltoc, "Very Low" = mbajoc),
                      labels = c("Very Low", "Low", "Medium", "High", "Very High"))+
  # coord_fixed()+
  # geom_abline()+
  # scale_x_continuous(name = "percentage language",limits = c(0,1))+
  # scale_y_continuous(name = "percentage mathematics",limits = c(0,100))+
  ggtitle("Achievement Level IV by marginalization")+
  th+
  theme(
    axis.text.x = element_text(angle = 0),
    panel.grid.minor = element_blank(),
    legend.position = "none")

mt + facet_grid(cols=vars(gm)) +
    theme(strip.background = element_rect(fill="#b53a31"))

# ggsave("~/Documents/UVM/courses/DS2/assignments/HW6/math_vs_language_IV.png",height = 20, width = 40, units = "cm")

Maps and things

I defined a function to get the geometries in an easy way. TODO: update fortify with broom

indig_pop <- indig %>% left_join(conapo,by="CLAVE")
indig_pop %<>% mutate(p_indig = IPOB_INDI/POP)
municipios_pop <- geom_muni %>% left_join(indig_pop, by = c("CVEGEO"="CLAVE"))
geom_muni$CVEGEO %>% unique() %>% length()

We have 7 municipalities without indigenous information, I don’t know if it is because they were instituted between 2015 and 2018, which is the last CONAPO available data and the geometry data that I have. (2456 to 2463, as in 2020 there are 2464)

ggplot()+
  geom_polygon(data = municipios_pop,
               aes(long,lat,label=CVEGEO,group=group,
                   fill=p_indig*100),color="grey",size=0.1)+
  coord_map(projection = "mercator")+
  scale_fill_continuous(name="Percentage",low="white", high="#0f5a5e",
                       guide="colorbar",limits=c(0,100))+
  labs(title="Indigenous Population")+
  theme(#panel.grid.major.y = element_line(color = "gray"),
        panel.grid.major.y = element_blank(),
        text = element_text(color = "gray20"),
        panel.background = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        legend.box = "horizontal",
        legend.text = element_text(size = 10),
        plot.caption = element_text(hjust=0),
        plot.title = element_text(size = 16, face = "bold",hjust=0.39))

ggsave("../../figs/map_indigenous.png",height = 30, width = 30, units = "cm")

municipios_nal <- geom_muni %>% left_join(mar_scores, by = c("CVEGEO"="CLAVE"))

A tilt in the map allows a better stacking

# Shear/scale matrix [[2,1],[0,1]] obtained by some trial and error:
sm <- matrix(c(1,-0.3,0,1),2,2)
# Get transformed coordinates:
xy <- as.matrix(municipios_nal[,c("long","lat")]) %*% sm
# Add xy as extra columns in fortified data:
municipios_nal$x <- xy[,1]; municipios_nal$y = xy[,2]

min(municipios_nal$x)
max(municipios_nal$x)

max(municipios_nal$y) - min(municipios_nal$y)
min(municipios_nal$y)
mean(municipios_nal$y)

Stacked map for the 4 levels of math

and language

p <- ggplot() +
  geom_polygon(data = municipios_nal, aes(x, y,label=CVEGEO, group=group,fill=(I_len)),color="white",size=0.01)+
  annotate("text",x=min(municipios_nal$x)+10,y=mean(municipios_nal$y),color="gray35",label="Level I")+
  geom_polygon(data = municipios_nal, aes(x, y-14,label=CVEGEO, group=group,fill=(II_len)),color="white",size=0.01)+
  annotate("text",x=min(municipios_nal$x)+10,y=mean(municipios_nal$y-14),color="gray35",label="Level II")+
  geom_polygon(data = municipios_nal, aes(x, y-28,label=CVEGEO, group=group,fill=(III_len)),color="white",size=0.01)+
  annotate("text",x=min(municipios_nal$x)+10,y=mean(municipios_nal$y-28),color="gray35",label="Level III")+
  geom_polygon(data = municipios_nal, aes(x, y-42,label=CVEGEO, group=group,fill=(IV_len)),color="white",size=0.01)+
  annotate("text",x=min(municipios_nal$x)+10,y=mean(municipios_nal$y-42),color="gray35",label="Level IV")+
  labs(title="Language")+
  scale_fill_continuous(name = "Mean percentage",low="white", high="#5a5a0d", guide="colorbar")+#,na.value="#f2f2f2")+
  coord_map(projection = "mercator")+
  theme(#panel.grid.major.y = element_line(color = "gray"),
        panel.grid.major.y = element_blank(),
        text = element_text(color = "gray20"),
        panel.background = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        # legend.box = "horizontal",
        legend.text = element_text(size = 10),
        plot.caption = element_text(hjust=0),
        plot.title = element_text(size = 16, face = "bold",hjust=0.39))

ggsave("../../figs/stacked_language.png",p,height = 30, width = 30, units = "cm")

I want to compute the density of schools by square kilometer

query <- load_table(con,public,geom_muni)
muni_df <- query %>% retrieve_result()
clear_results(con)

Build the polygons and compute the area in Km^2

polman <- lapply(muni_df$WKT,rgeos::readWKT)
rgeos::gUnion(polman[[1]],polman[[2]])
length(polman)

the_areas <- sapply(polman,geosphere::areaPolygon)
areas_df <- tibble(CVEGEO=muni_df$CVEGEO,area=the_areas/1000000)

Now we need the number of Schools but not from the test, we’ll use the data from Consejo Nacional de Fomento Educativo (CONAFE) through the Censo de Escuelas, Maestros y Alumnos de Educación Básica y Especial (CEMABE) Nop, now I now that CONAFE is a program to contribute children in highly marginalized areas to conclude their basic comunitary education, with economical support to educational figures and scholarly supplies to students. It would be interesting to analize this into deep!!!

This is the number of schools in the planea set

unique(planea$clave) %>% length()

schools_mun <- schools %>% group_by(CLAVE) %>% summarise(num=n())
areas_df %<>% left_join(schools_mun,by=c("CVEGEO"="CLAVE"))
areas_df %<>% mutate(density = num/area)

The density is log-normal As well as the test results. Should consider this when doing regressions.

ggplot(areas_df)+
  geom_histogram(aes(log10(density)))+
  th

geom_mun_school <- geom_muni %>% left_join(areas_df, by = c("CVEGEO"))

library(latex2exp)
ggplot()+
  geom_polygon(data = geom_mun_school,
               aes(long,lat,label=CVEGEO,group=group,
                   fill=(density)),color="grey",size=0.01)+
  coord_map(projection = "mercator")+
  # scale_fill_continuous(name="Percentage",low="white", high="#4c1262", guide="colorbar")+#,limits=c(0,100))+
  scale_fill_gradient(trans="log10",name="Density",low="white", high="#4c1262", guide="colorbar",labels=c(0,0.5,1,1.5,2,2.5))+#,limits=c(0,100))+
  labs(title=TeX("Schools by $Km^2$"))+
  theme(#panel.grid.major.y = element_line(color = "gray"),
        panel.grid.major.y = element_blank(),
        text = element_text(color = "gray20"),
        panel.background = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        legend.box = "horizontal",
        legend.text = element_text(size = 10),
        plot.caption = element_text(hjust=0),
        plot.title = element_text(size = 16, face = "bold",hjust=0.39))

ggsave("../../figs/map_school_density.png",height = 30, width = 30, units = "cm")

I think I haven’t written the ingestion code for the schools into the sql database. TODO Although I’m not using it.

conafe <- read_delim("/home/ollin/Documents/UVM/courses/DS2/project/planea/data/clean/tr_conafe.csv",delim="|")
conafe %<>% mutate(CVEGEO=paste0(ENT,MUN))
conafe %>% group_by(CVEGEO) %>% summarise(num = n())

Maps are cool but we need other stuff, let’s make ridges for achievement

library(ggridges)

mar_tidy <- mar_scores %>% pivot_longer(c(I_len,II_len,III_len,IV_len,I_mat,II_mat,III_mat,IV_mat)) %>% unique()

First for all

ggplot(mar_tidy[grep("len",mar_tidy$name),], aes(x = value, y = name, height = stat(density))) +
  geom_density_ridges(stat = "density")+
  scale_x_continuous(name="Percentage")+
  scale_y_discrete(name="Level")+
  th
ggsave("../../figs/ridges_len.png",height = 30, width = 30, units = "cm")

ggplot(mar_tidy[grep("mat",mar_tidy$name),], aes(x = value, y = name, height = stat(density))) +
  geom_density_ridges(stat = "density")+
  scale_x_continuous(name="Percentage")+
  scale_y_discrete(name="Level")+
  th
ggsave("../../figs/ridges_math_log.png",height = 30, width = 30, units = "cm")

We conducted an Anderson-Darling test to assess the normality of distributions along achievement levels

compute how likely are them to be drawned by a normal distribution

This means that our values are not normal

compute how likely are them to be drawned by a the same distribution

combos <- combn(4,2)
plyr::adply(combos, 2, function(x) {
  test <- t.test(Datas[, as.numeric(x[1])], Data[, as.numeric(x[2])],alternative = "t")

  out <- data.frame("var1" = colnames(Data)[x[1]]
                    , "var2" = colnames(Data[x[2]])
                    , "t.value" = sprintf("%.5f", test$statistic)
                    ,  "df"= test$parameter
                    ,  "p.value" = test$p.value#sprintf("%.5f", test$p.value)
                    )
  return(out)

})
# mar_scores %>%
#   dplyr::select(I_len,II_len,III_len,IV_len) %>%
#   dplyr::select_if(is.numeric) %>%
#   purrr::map_df(~ broom::tidy(t.test(. ~ grp)), .id = 'var')

ggplot(mar_tidy[grep("len",mar_tidy$name),], aes(x = value, y = name, height = stat(density))) +
  geom_density_ridges(stat = "density")+
  scale_x_log10(name="Percentage")+
  scale_y_discrete(name="Level")+
  th

ggsave("../../figs/ridges_language_log.png",height = 30, width = 30, units = "cm")

compute how likely are them to be have the same mean

compute how likely are them to have the same mean

combos <- combn(4,2)
plyr::adply(combos, 2, function(x) {
  test <- t.test(Data[, as.numeric(x[1])], Data[, as.numeric(x[2])],alternative = "t")

  out <- data.frame("var1" = colnames(Data)[x[1]]
                    , "var2" = colnames(Data[x[2]])
                    , "t.value" = sprintf("%.5f", test$statistic)
                    ,  "df"= test$parameter
                    ,  "p.value" = test$p.value#sprintf("%.5f", test$p.value)
                    )
  return(out)

})
# mar_scores %>%
#   dplyr::select(I_len,II_len,III_len,IV_len) %>%
#   dplyr::select_if(is.numeric) %>%
#   purrr::map_df(~ broom::tidy(t.test(. ~ grp)), .id = 'var')

Then the facets by marginalization level

ggplot(mar_tidy[grep("mat",mar_tidy$name),], aes(x = value, y = name,fill=gm,alpha=0.6, height = stat(density))) +
  scale_fill_manual(name = "Marginalization", values = c("High" = altoc, "Low" = bajoc, "Medium" = medioc, "Very High" = maltoc, "Very Low" = mbajoc),
                      labels = c("Very Low", "Low", "Medium", "High", "Very High"))+
  geom_density_ridges(stat = "density")+
  scale_x_log10()+
  th+
  theme(
    axis.text.x = element_text(angle = 0),
    panel.grid.minor = element_blank(),
    legend.position = "none")+
  facet_grid(~gm)+
  theme(strip.background = element_rect(fill="#b53a31"))

ggplot(mar_tidy[grep("len",mar_tidy$name),], aes(x = value, y = gm,fill=name,alpha=0.6, height = stat(density))) +
  scale_fill_manual(name = "Marginalization", values = c("High" = altoc, "Low" = bajoc, "Medium" = medioc, "Very High" = maltoc, "Very Low" = mbajoc),
                        labels = c("Very Low", "Low", "Medium", "High", "Very High"))+
    geom_density_ridges(stat = "density")+
  scale_x_log10()+
    th+
    theme(
      axis.text.x = element_text(angle = 0),
      panel.grid.minor = element_blank(),
      legend.position = "none")+
    facet_grid(~name)+
    theme(strip.background = element_rect(fill="#b53a31"))

ggplot(mar_tidy[!grepl("II_",mar_tidy$name) & grepl("I_",mar_tidy$name),], aes(x = value, y = name,fill=gm,alpha=0.6, height = stat(density))) +
  scale_fill_manual(name = "Marginalization", values = c("High" = altoc, "Low" = bajoc, "Medium" = medioc, "Very High" = maltoc, "Very Low" = mbajoc),
                        labels = c("Very Low", "Low", "Medium", "High", "Very High"))+
    geom_density_ridges(stat = "density")+
  scale_x_log10(name="Percentage")+
  scale_y_discrete(name=" ")+
    th+
    theme(
      axis.text.x = element_text(angle = 0),
      panel.grid.minor = element_blank(),
      legend.position = "none")+
    facet_grid(~gm)+
    theme(strip.background = element_rect(fill="#b53a31"))

ggsave("../../figs/ridges_marg_I.png",height = 15, width = 30, units = "cm")

We test the homogeneity of variances with an F-test because we want to see how well does an achievement level in one test explains the other.

tocomp <- mar_scores %>%
  group_by(gm) %>%
  select(gm,I_len,I_mat,IV_len,IV_mat)
logss <- sapply(tocomp %>% ungroup(gm) %>% select(-gm), function(x) log10(as.numeric(x)) )
tocomp[2:5] <- logss
tocomp <- tocomp[apply(tocomp %>% ungroup(gm) %>% select(-gm), 1, function(row) all(is.finite(row))),]
# tocomp <- tocomp %>% sample_n(250,replace = F)

tocomp %>%
  group_by(gm) %>%
  summarise(pval = var.test(I_len,I_mat,alternative="t")$p.value)

ggplot(mar_tidy[grepl("IV_",mar_tidy$name),], aes(x = (value), y = name,fill=gm,alpha=0.6, height = stat(density))) +
  scale_fill_manual(name = "Marginalization", values = c("High" = altoc, "Low" = bajoc, "Medium" = medioc, "Very High" = maltoc, "Very Low" = mbajoc),
                        labels = c("Very Low", "Low", "Medium", "High", "Very High"))+
    geom_density_ridges(stat = "density")+
  scale_x_log10(name="Percentage")+
  scale_y_discrete(name=" ")+
    th+
    theme(
      axis.text.x = element_text(angle = 0),
      panel.grid.minor = element_blank(),
      legend.position = "none")+
    facet_grid(~gm)+
    theme(strip.background = element_rect(fill="#b53a31"))
ggsave("../../figs/ridges_marg_IV.png",height = 15, width = 30, units = "cm")

Use the Fligner-Killen median test to compare the homogeneity of variances

tocomp <- mar_scores %>%
  group_by(gm) %>%
  select(gm,I_len,I_mat,IV_len,IV_mat)
logss <- sapply(tocomp %>% ungroup(gm) %>% select(-gm), function(x) log10(as.numeric(x)) )
tocomp[2:5] <- logss
tocomp <- tocomp[apply(tocomp %>% ungroup(gm) %>% select(-gm), 1, function(row) all(is.finite(row))),]
# tocomp <- tocomp %>% sample_n(300,replace = F)

tocomp %>%
  group_by(gm) %>%
  summarise(pval = var.test(IV_len,IV_mat,alternative="t")$p.value)

tocomp <- schools[!is.na(schools$gm),] %>%
  group_by(gm) %>%
  select(gm,I_porc_len,I_porc_mat,IV_porc_len,IV_porc_mat)
logss <- sapply(tocomp %>% ungroup(gm) %>% select(-gm), function(x) log10(as.numeric(x)) )
tocomp[2:5] <- logss
tocomp <- tocomp[apply(tocomp %>% ungroup(gm) %>% select(-gm), 1, function(row) all(is.finite(row))),]
# tocomp <- tocomp %>% sample_n(250,replace = F)

tocomp %>%
  group_by(gm) %>%
  summarise(pval = var.test(IV_porc_len,IV_porc_mat,alternative="t")$p.value)

But what about school type

mar_tidy <- mar_scores %>% pivot_longer(c(I_len,II_len,III_len,IV_len,I_mat,II_mat,III_mat,IV_mat)) %>% unique()

mar_scores %<>% left_join(schools %>% select(CLAVE,tipo_escuela),by="CLAVE") %>% unique()
mar_tidy2 <- mar_scores %>% pivot_longer(c(I_len,II_len,III_len,IV_len,I_mat,II_mat,III_mat,IV_mat)) %>% unique()
ggplot(mar_tidy2[grepl("I_",mar_tidy2$name),], aes(x = (value), y = name,fill=tipo_escuela,alpha=0.6, height = stat(density))) +
  # scale_fill_manual(name = "Marginalization", values = c("High" = altoc, "Low" = bajoc, "Medium" = medioc, "Very High" = maltoc, "Very Low" = mbajoc),
  #                       labels = c("Very Low", "Low", "Medium", "High", "Very High"))+
    geom_density_ridges(stat = "density")+
  scale_x_log10()+
    th+
    theme(
      axis.text.x = element_text(angle = 0),
      panel.grid.minor = element_blank(),
      legend.position = "none")+
    facet_grid(~tipo_escuela)+
    theme(strip.background = element_rect(fill="#b53a31"))

Build the distribution difference variable

schools$I_porc_mat - schools$I_porc_len
schools %<>% mutate(diffI = I_porc_mat - I_porc_len, diffIV = IV_porc_mat - IV_porc_len)
sch_tidy <- schools %>% pivot_longer(c(diffI, diffIV)) %>% unique()

ggplot(sch_tidy, aes(x = (value), y = name,fill=gm,alpha=0.6, height = stat(density))) +
  scale_fill_manual(name = "Marginalization", values = c("High" = altoc, "Low" = bajoc, "Medium" = medioc, "Very High" = maltoc, "Very Low" = mbajoc),
                        labels = c("Very Low", "Low", "Medium", "High", "Very High"))+
    geom_density_ridges(stat = "density")+
  scale_x_log10(name="Percentage")+
  scale_y_discrete(name=" ")+
    th+
    theme(
      axis.text.x = element_text(angle = 0),
      panel.grid.minor = element_blank(),
      legend.position = "none")+
    facet_grid(~gm)+
    theme(strip.background = element_rect(fill="#b53a31"))

ggplot(schools[complete.cases(schools$gm),],aes(diffIV,fill=gm))+
  geom_density()

ggplot(schools[!is.na(schools$gm),],aes(x=(diffIV+100)/200,y=gm,fill=gm,height = stat(density)))+
  geom_density_ridges(stat="density")+
  scale_x_log10()


# ggsave("../../figs/ridges_marg_IV.png",height = 15, width = 30, units = "cm")

ggplot(schools[!is.na(schools$gm),],aes(x=(diffI+100)/200,y=gm,fill=gm,height = stat(density)))+
  geom_density_ridges(stat="density")+
  scale_x_log10()#limits = c(-50,50))

mate - lengua

ci_median <- function(x){
  bootmed = apply(matrix(sample(x, rep=TRUE, 10^4*length(x)), nrow=10^3), 1, median)
  quantile(bootmed, c(.025,0.5, 0.975))
}

schools[complete.cases(schools$gm),] %>% group_by(gm) %>% summarise(median(diffI),median(diffIV))
schools %>% group_by(gm) %>% summarise(ci_median(log10((diffIV+100)/200))[2])

schools[complete.cases(schools$gm),] %>% group_by(gm) %>% summarise(median(diffI),ci_median(diffIV)[1],ci_median(diffIV)[3])

Let’s check the Bayesian networks

Puerto Morelos became municipality in 2011

encode_ordinal <- function(x, order = unique(x)) {
  x <- as.numeric(factor(x, levels = order, exclude = NULL))
  x
}

schools <- schools %>% select(-ovsd,-ovsdse,-ind0a100)
schools <- schools[complete.cases(schools),]
areas_df <- areas_df %>% mutate(CLAVE=CVEGEO) %>% select(-CVEGEO)
indig

for_bn <- schools %>% left_join(indig %>% select(CLAVE,IPOB_INDI,IPhli))
for_bn <- for_bn %>% mutate(p_indig = IPOB_INDI/pop_marg,lang_indig=IPhli/IPOB_INDI)
for_bn %<>% left_join(areas_df,by="CLAVE")

for_bn[!complete.cases(for_bn),]$IPOB_INDI <- 0
for_bn[!complete.cases(for_bn),]$p_indig <- 0
databn <- for_bn %>% ungroup(CLAVE,NOM_ENT) %>%
  select(diffI,diffIV,turno,gm,tipo_escuela,illiteracy,elementary,no_sewage,no_electricity,no_water,overcrowding,dirt_floor,less_5k,less_2minwage,im,lugar,lugarest,IPOB_INDI,p_indig,area,num,density,lang_indig)
  # select(I_porc_len,II_porc_len,III_porc_len,IV_porc_len,I_porc_mat,II_porc_mat,III_porc_mat,IV_porc_mat,turno,gm,tipo_escuela,illiteracy,elementary,no_sewage,no_electricity,no_water,overcrowding,dirt_floor,less_5k,less_2minwage,im,lugar,lugarest,IPOB_INDI,p_indig,area,num,density,lang_indig)

databn

# I should automate this
databn$turno <- encode_ordinal(databn$turno)
databn$gm <- encode_ordinal(databn$gm)
databn$tipo_escuela <- encode_ordinal(databn$tipo_escuela)

# databn$turno <- as.double(databn$turno)
# databn$gm <- as.double(databn$gm)
# databn$num <- as.double(databn$num)
# databn$tipo_escuela <- as.double(databn$tipo_escuela)

databn2 <- databn

databn$I_porc_len <- log10(databn$I_porc_len)
databn$II_porc_len <- log10(databn$II_porc_len)
databn$III_porc_len <- log10(databn$III_porc_len)
databn$IV_porc_len <- log10(databn$IV_porc_len)

databn$I_porc_mat <- log10(databn$I_porc_mat)
databn$II_porc_mat <- log10(databn$II_porc_mat)
databn$III_porc_mat <- log10(databn$III_porc_mat)
databn$IV_porc_mat <- log10(databn$IV_porc_mat)

databn$density <- log10(databn$density)
databn$p_indig <- log10(databn$p_indig)
databn$lang_indig <- log10(databn$lang_indig)

databn_finite <- databn[apply(databn, 1, function(row) all(is.finite(row))),]
databn_finite <- databn[apply(databn %>% select(-turno,-gm,-num,-tipo_escuela), 1, function(row) all(is.finite(row))),]
databn_finite %<>% select(-im,-lugar,-lugarest,-IPOB_INDI,-area,-num)

Delete the highly correlated variables

set.seed(42)
# corrM <- cor(databn_finite[,c(3:17)])
corrM <- cor(databn_finite)
highlyCorrelated <- findCorrelation(corrM, cutoff=0.70)
colnames(databn_finite[,c(3:17)])[highlyCorrelated]
databn_finite %<>% select(-illiteracy,-elementary,-overcrowding,-dirt_floor,-less_5k)
databn_finite %<>% select(-gm)

Variable selection

# control <- trainControl(method="cv", number=10)
# model <- caret::train(I_porc_len~., data=databn_finite, method="cforest", preProcess="scale", trControl=control)
# # estimate variable importance
# importance <- varImp(model, scale=FALSE)
# # summarize importance
# print(importance)
# # plot importance
# plot(importance)

# set.seed(42)
# databn_finite2 <- as.data.frame(databn_finite)
# #Controlling the feature selection, using default functions, with 10 folds cross validation
# control <- rfeControl(functions=rfFuncs, method="cv", number=10)
# # run the RFE algorithm
# # using recursive feature elimination or backwards selection
# # let's check first if it works the same for each level and test
# results <- rfe(x=databn_finite2[,c(1,10:19)], y=databn_finite2[,2], rfeControl=control)
# # results <- rfe(databn[,c(1,10:27)], databn[,n], sizes=c(1:18), rfeControl=control)
# # summarize the results
# print(results)
# # list the chosen features
# predictors(results)
# # plot the results
# plot(results, type=c("g", "o"))

I haven’t figured how to save this graphs with the chunk output inline in a notebook, so we have to dissable it to run and save the plots

for (i in 1:7){
  bnf <- bnlearn::h2pc(databn_finite[,c(i,9:17)])
  png(paste0("/figs/graph_",i,".png"))
  plot(graphviz.plot(bnf, shape = "ellipse"))
  dev.off()
  str.diff <- bnlearn::boot.strength(databn_finite[,c(i,9:17)],R=20,algorithm="h2pc")
  avg.diff <- averaged.network(str.diff)#,threshold = 0.7)
  png(paste0("/figs/graph_",i,"_boots.png"))
  plot(graphviz.plot(avg.diff, shape = "ellipse"))
  dev.off()
}

# the_columns <- c(diffI, turno, tipo_escuela, illiteracy, elementary, no_sewage, no_electricity, no_water, overcrowding, dirt_floor, less_5k, less_2minwage, im, p_indig, density, lang_indig)
# bnf <- bnlearn::h2pc(databn_finite[,c(1,3:11)])
bnf <- bnlearn::h2pc(databn_finite %>% select(diffI, turno, tipo_escuela, illiteracy, elementary, no_sewage, no_electricity, no_water, overcrowding, dirt_floor, less_5k, less_2minwage, im, p_indig, density, lang_indig))

# str.diff <- bnlearn::boot.strength(databn_finite[,c(1,3:11)],R=20,algorithm="h2pc")
str.diff <- bnlearn::boot.strength(databn_finite %>% select(diffI, turno, tipo_escuela, illiteracy, elementary, no_sewage, no_electricity, no_water, overcrowding, dirt_floor, less_5k, less_2minwage, im, p_indig, density, lang_indig),R=20,algorithm="h2pc")

avg.diff <- averaged.network(str.diff)#,threshold = 0.7)
graphviz.plot(avg.diff, shape = "ellipse")

bnf <- bnlearn::h2pc(databn_finite[,c(2,3:11)])
  str.diff <- bnlearn::boot.strength(databn_finite[,c(2,3:11)],R=20,algorithm="h2pc")
  avg.diff <- averaged.network(str.diff)#,threshold = 0.7)
  graphviz.plot(avg.diff, shape = "ellipse")

bnf <- bnlearn::h2pc(databn_finite[,c(1,9:17)])
fitted <- bn.fit(bnf,databn_finite[,c(1,9:17)])

graphviz.plot(bnf, shape = "ellipse")

bnf <- bnlearn::h2pc(databn_finite[,c(1,3:11)])
bnf <- bnlearn::reverse.arc(bnf,from="tipo_escuela",to="IV_porc_len")
bnf <- bnlearn::reverse.arc(bnf,from="turno",to="IV_porc_len")
bnf <- bnlearn::reverse.arc(bnf,from="density",to="IV_porc_len")
fitted <- bn.fit(bnf,databn_finite[,c(4,3:11)])

graphviz.plot(bnf, shape = "ellipse")

bnf <- bnlearn::h2pc(databn_finite[,c(2,3:11)])
bnf <- bnlearn::reverse.arc(bnf,from="tipo_escuela",to="I_porc_mat")
bnf <- bnlearn::reverse.arc(bnf,from="turno",to="I_porc_mat")
bnf <- bnlearn::reverse.arc(bnf,from="overcrowding",to="I_porc_mat")
bnf <- bnlearn::reverse.arc(bnf,from="lang_indig",to="I_porc_mat")
fitted <- bn.fit(bnf,databn_finite[,c(2,3:11)])

graphviz.plot(bnf, shape = "ellipse")

bnf <- bnlearn::h2pc(databn_finite[,c(8,9:17)])
bnf <- bnlearn::reverse.arc(bnf,from="tipo_escuela",to="IV_porc_mat")
bnf <- bnlearn::reverse.arc(bnf,from="turno",to="IV_porc_mat")
fitted <- bn.fit(bnf,databn_finite[,c(8,9:17)])

bn.fit.qqplot(fitted)
bn.fit.barchart(fitted$IV_porc_mat)
bn.cv(databn_finite[,c(1,9:17)],bn="h2pc")
bn.cv(databn_finite[,c(4,9:17)],bn="h2pc")

bn.cv(databn_finite[,c(5,9:17)],bn="h2pc")
a<-bn.cv(databn_finite[,c(8,9:17)],bn="h2pc")

graphviz.plot(bnf, shape = "ellipse")

1 - mat log likelihood 14.78539
4 - mat log likelihood 15.17735

1 - lan log likelihood 14.93725
4 - lan log likelihood 14.79238